Plane waves in periodic, quadratically nonlinear slab waveguides: stability and exact 

Fourier structure 
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We consider the propagation of broad optical beams through slab waveguides with a purely 
quadratic nonlinearity and containing linear and nonlinear long-period quasi-phase-matching grat- 
ings. An exact Floquet analysis on the periodic, plane-wave solution shows that the periodicity 
can drastically alter the growth rate of the modulational instability but that it never completely 
removes the instability. The results are confirmed by direct numerical simulation, as well as through 
a simpler, approximate theory for the averaged fields that accurately predicts the low-frequency part 
of the gain spectrum. OCIS codes: 190.5530, 230.4320, 190.4410, 190.5940, 



I. INTRODUCTION 

Optical materials with \^ nonlinearity offer many ad- 
vantages as candidates for all-optical processing applica- 
tions. For example, they offer a fast and strong nonlin- 
ear response, and through cascading \^ : x' 2 ' processes 
, they can support stable localized solitary- wave solu- 
tions in more than one dimension J2|-^| , in contrast with 
systems that have only a cubic nonlinear response. Ex- 
perimental studies have confirmed the existence of one- 
[7\ and two-dimensional Q spatial and spatio-temporal 
[3| quadratic solitons and have demonstrated that t hey 
act as dynamically reconfigurable guiding structures [10| , 
with switching behavior controllable through phase- 
and intensity-dependent |L2) effects. Even though special 
'dark-vortex' solitons have been generated in quadratic 
media jljj, coherent dark solitons are generally unsta- 
ble in homogeneous materials with only quadratic non- 
linearity, due to a modulational instability (MI) in the 
background plane wave. 

Like solitons, MI occurs due to an interplay between 
nonlinear and dispersive or diffractive effects. In non- 
linear optics, it appears as a transverse instability that 
breaks up a broad optical beam [jl4} , thereby acting as a 
precursor for the formation of stable bright spatial soli- 
tons, as has been confirmed experimentally p5| . Con- 
versely, the stable propagation of dark solitons relies on 
the stability of the constant-intensity background and 
thus requires the absence of MI. The mismatch Q and in- 
tensity |16| dependence of the instability period has been 
investigated experimentally, as has the intensity and fre- 
quency dependence of the gain p7| . 

In all materials, there is of course an inherent cu- 
bic (Kerr) nonlinearity. The combined nonlinearities lead 
to novel effects, such as the existence of bright solitons 
in an intrinsically defocusing cubic medium |lqj . Also, if 
the defocusing cubic nonlinearity is sufficiently strong, it 
may eliminate MI ]19fl . Unfortunately, the cubic nonlin- 
earity in conventional \^ materials is usually focusing 
and thus of the wrong type for these effects to occur. 

However a cubic nonlinearity of a different kind may 
be induced through long-period variations in the linear 



or nonlinear susceptibilities. Such periodicities are used 
for quasi-phase-matching (QPM), which compensates the 
material mismatch between the fundamental and second- 
harmonic (SH) modes. The lowest-order effect of QPM 
is to induce an effectively phase-matched quadratic non- 
linearity in the equations for the averaged fields. To first 
order, the non-phase-matched oscillatory components of 
the fields induce cubic nonlinearities in the 

same manner as the non-phase-matched interactions do 
in the cascading limit |24j| . For certain materials and 
etching techniques, altering the nonlinear susceptibility 
also induces a change in the linear susceptibility. The 
action of such a simultaneous linear grating can be to re- 
duce the effective quadratic nonlinearity [p5|-p9t , making 
the effective nonlinearity more cubic, as in the cascading 
limit. 

In this paper, we investigate the effects of such peri- 
odicities on the MI properties of beams in x^ materi- 
als, considerably extending the initial results presented in 
Ref. |p3. We provide a comprehensive and comparative 
treatment of the different approaches (exact and approx- 
imate) and use a variety ofparticular cases to illustrate 
general properties. Section y introduces the system stud- 
ied in the res t of the paper. In analyzing this system, 
we first (Sec. Ill) seek modulationally stable dark soli- 
tons in an approximate averaged system. Previous work 
p2| , p3| has confirmed the accuracy of the averaged theory 
in describing the average properties of solitons and CW 
switching J3(j|nj , showing that this intrinsically quadrat- 
ically nonlinear system behaves as though it were gov- 
erned by cubic as well as quadratic nonlinearities. The 
theory also describes space-varying behavior, but to a 
lower accuracy. We discover that the cubic terms that 
result from particular gratings can be large enough and 
of the correct sign to prevent MI. 

The second part of this work compares the stability 
predictions of the averaged theory with those of exact cal- 
culations. The exact methods involve finding the exact 
periodic plane-wave solutions in the full (nonaveraged) 
system, with the higher harmonics included to all orders 
(Sec. IV). We apply a Floquet analysis to these periodic 
solutions to determine the growth rate of linear pertur- 
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bations (Sec. 
experiments | 



which can now be directly probed in 
Additionally, we directly simulate the 
periodic-field equations (Sec. VI), using for initial condi- 
tions the averaged solutions seeded with noise. 

According to the exact calculations, the periodicities 
can drastically alter the MI gain spectrum, but they 
do not entirely remove the inherent instability. In the 
regime of efficient QPM, we find that the MI gain spec- 
trum contains two distinct and well-separated features 
with fundamentally different physical origins (see Sec. 
VB). The low- frequency part of the gain is accurately 



predicted by the averaged theory and disappears for cer- 
tain grating modulations. The high-frequency part of the 
spectrum is related to the inherent gain of the non-phasc- 
matchcd material (i.e., with no gratings) and appears to 
be unavoidable. However, because they are consistently 
small, the high-frequency peaks can be ignored under a 
less stringent definition of experimental stability (see Sec. 
VP ). This gives the possibility of stable plane waves as 
predicted by the averaged theory. 



II. SYSTEM 

We consider a CW beam (carrier frequency lo) inter- 
acting with its second harmonic (SH), both propagat- 
ing in a lossless ID slab waveguide under conditions 
for type I second-harmonic generation (SHG). Both the 
linear and quadratic nonlinear susceptibilities are peri- 
odic along the Z-direction of propagation but are con- 
stant along the transverse X-direction. We assume a 
weak modulation of the refractive index (Anj(Z)/fij<^l, 



where nj(Z) 



- Arij(Z) and j refers to the frequency 



juj). We consider gratings for forward QPM only. The 
grating period is then much longer than the optical pe- 
riod, in which case Bragg reflections can be neglected. 
The evolution of the slowly varying beam envelopes is 
then described by Q 

id z E 1 + \d 2 x E x + ai Ex + xE{E 2 exp (if3Z) = 0, 

id z E 2 + \d 2 x E 2 + 2a 2 E 2 + X E\ exp {-i/3Z) = 0, (1) 

where Ei=Ei(X, Z) and E 2 =E 2 (X, Z) are the envelope 
functions of the fundamental and SH, respectively. The 
coordinates X and Z are in units of the input beam 
width Xq and the diffraction length Ld=k\X^ respec- 
tively. The parameter (3—AkLd is proportional to the 
mismatch Ak=k 2 — 2ki, where kj=jwhj/c is the average 
wave number. Thus (3 is positive for normal dispersion 
and negative for anomalous dispersion. The normalized 
refractive- index grating is aj(Z)=LdwAnj(Z)/c, and the 
normalized nonlinear grating is x(Z) = Ed UJ d c g (Z) / (fiic) , 
where d c g /2 is given in MKS units. The model 
(|l|) describes both temporal and spatial solitons p3| , ^4[ . 
Throughout this paper, we focus on first-order QPM us- 
ing conventional square gratings with a 50% duty cycle, 
as shown in Fig. [|. 
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FIG. 1. Periodic linear and quadratically nonlinear grat- 
ings, a.j(Z) and x(Z), with period 2Lo=2n /\k g \. 

To consider plane wave solutions with longitudi- 
nal wave number offset A and transverse wave num- 
ber f2, we transform to the new variables w(x,z) = 
E 1 (X,Z)exp(iQ), v(x,z) = E 2 (X, Z) exp(2i6 + ifiZ), 
x = ^/rj(X + QZ), and z = r/Z, where Q(z) = Q.X - 
AZ - j at(z)dz and where n = |A + \Vt 2 \ for A ^ -\Vt 2 
and rj=\ otherwise. We have also introduced a residual 
effective mismatch f3 = (3 — k g , where k g is the spatial 
frequency of the gratings. This gives the rescaled equa- 
tions 



id z w + 2®x w — rw + x' w * v ex P {inz) = 
id z v + -^d 2 v — av + av + xw 2 exp (— inz) = 0, (2) 



where r = sgn(A + \Q 2 ), k = k g /rj, and a = 2r — f3' , 
for 0' = 0/r). Note that r=0 when A=-±fl 2 . This scal- 
ing by rj may be absorbed into the original normalization 
[Eq. (||)] through a redefinition of X (and hence L4) to 
be something other than the true beam width (diffrac- 
tion length). We choose to retain 77 explicitly, however, 
to show the connection with the standard normalization. 

To find solutions in this periodic system, we Fourier- 
expand the rescaled gratings, a' {z)=2(a 2 (Z) — a\{Z))/rj 
and x'{z) = x{Z)/w- 

a = a g n exp (innz), x = ^0 + df g n exp (innz), 



(3) 



where a' = 2{a 2 — a\)/rj and where g n — 2s/ {inn) for 
n odd and g n = for n even. We include s = sgn(fc) 
in the Fourier expansions to ensure that the expansions 
give the same physical grating form for both signs of k. 
Plane waves in system ([l]) and (^) have an intensity that 
is constant in the transverse direction but periodic in Z, 
driven by the periodic variation of the susceptibilities. 
Thus they have the same spatial frequency k as the grat- 
ings, and we expand them in Fourier series also: 
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w = w n{ x 3 z ) ex P {intzz), v — v n (x, z) exp (intzz). 

n n 

(4) 

An accurate treatment of the problem requires that all 
Fourier components (up to some cutoff) be included. 
This we do in Sec. EyL However we first present an 
approximate theory for the DC components (w, v) — 
(wq,vq), which suffices for many calculations. 

III. AVERAGED-FIELD THEORY 

A. Plane-wave solutions 

Three physical length scales are in play: the diffrac- 
tion length Ld, the coherence length L c , and the grat- 
ing domain length L . In our rescaled units = f], 
L c = ??7r/|/3| = 7r/|/3'|, and L = tt/\k\. We assume a 
typical QPM grating with a domain length that is much 
shorter than the diffraction length and use perturbation 
theory ^] with the small parameter e = Lo/Ld <C 1 . We 
insert the Fourier expansions (||) and (|j) into the dynam- 
ical equations and assume the harmonics w n ^o and v n ^Q 
to be of order e. Furthermore, efficient phase matching 
is assumed, with the domain length being close to the 
coherence length, L ~ L c , so the residual mismatch is 
small, \(3'\ <C \k\- To first order (e 1 ), this gives the har- 
monics 

(d'g n -i + d' 6 n ,i)wQV 
w n ^a = , 

UK 

a'g n v + (d'g n+1 + d' 5 n -i)w% 

Vn^O = ■ , (5) 

nn 

where 8 n ^ m is Kronecker's delta and where we have as- 
sumed that the coefficients w n (x, z) and v n (x,z) vary 
slowly in z compared to exp(mz). Using these solutions, 
we obtain to first order the averaged equations for the 
DC components (w,v) = (wq,Vo)- 

id z w + 2®%™ ~ r w + pw*v + 7(|w| 2 — |u>| 2 )u) = 0, 
id z v + -^d^v — ay + p*w 2 + 2j\w\ 2 v = 0, (6) 

These equations also describe mth order QPM (where 
P — [3 — mkg is ideally zero) and any other type of long- 
period grating. Incorporating time or the spatial y co- 
ordinate is also straightforward. The effective quadratic 
(p) and induced cubic (7) nonlinearities are simply given 
by sums over the Fourier coefficients |2(]], and for the 
square grating (^) are p = i2s {d' — d' a a' / n} /it and 
7= K 2 +d' 2 (f-8/7r 2 )]/ K . 

From Eqs. (ph follows the important result that cubic 
nonlinearities are induced by the nonlinear QPM grat- 
ings. This cubic nonlinearity has the form of asymmetric 
self- and cross-phase modulations (SPM and XPM) and 



is a result of non-phase-matched coupling between the 
wave at the main spatial frequency n and its higher har- 
monics. It is thus of a fundamentally different nature 
than the material Kerr nonlinearity, a fact which is re- 
flected in the the SPM being absent for the SH. 

The induced cubic nonlinearity depends only on the 
nonlinear grating and may lead to either a focusing or a 
defocusing effect, depending on the relative intensity of 
the fields and the sign of the phase mismatch (3' , since 
sgn(ft) — sgn (/?'). The effective nonlinearity has a 
contribution from the linear grating, which can either re- 
inforce or undermine the contribution from the nonlinear 
grating, depending on the physical situation [ p8p7| ]. 

Any plane wave will correspond to a stationary solu- 
tion (w S: v s ) of Eq. (](]), obtained by setting the deriva- 
tives to zero. By using a gauge transformation, we make 
w s real and positive. The solutions are then w s — y/y/\p\, 
v s = y/(<Jp — 2jyp), where y is a positive, real number 
that satisfies 

= 47 3 y 3 +7(1 + 47(r - a)) y 2 + a (ja - Ar^f - 1) y + ra 2 . 

(7) 

For a given r, 7 and er, there can thus be anywhere 
up to three solutions. The relative strength of the cu- 
bic nonlinearity is given by the parameter 7 = t/ | /O | 2 , 
which depends only on the material grating parameters 
k, a'/ (3', and d' /d'. Note that the residual-mismatch 
term a equals 2r for perfect phase matching. 
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FIG. 2. Ratio of intensities L = log (vf /wig) for plane-wave 
solutions, as a function of 7 and a, for r = — 1 

The solutions in parameter space (7, a) are illustrated 
by the intensity ratio R = v 2 /w 2 in Fig. || for r = — 1 . 
The figure shows that there are four extended branches 
of solution, which exist in three quadrants, the fourth 
quadrant being devoid of any solution. We label the per- 
vasive branch (G) and the upper branch in the negative 
quadrant (U2). The boundaries of (G) and (U2) lie along 
the 7 and a axes and are where the ratio R tends to in- 
finity The positive quadrant contains additional paired 
upper (Ul) and lower (LI) branches, which exist only 
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away from the axes, their boundary given by: 



B. Modulational instability 



cr(7) 



5 + 167 + 3(, 



(8) 



where f± = -1 + 80-y + I287 2 ± 16^/j(-l + 4^) 3 . Note 
that for 7 < 1/4, we must take the roots of f± that give 
the largest real a. Equation (||) tells us that (Ul) and 
(LI) exist only for 7 > and a > 2. 

Note that the r = — 1 solutions can be transformed 
to the r = 1 solutions by the replacement v s — * — V s , 
a — > —a and 7 — > —7. However, the MI properties differ. 
In addition to the solutions in Fig. [| there is also a triv- 
ial solution (w s ,v s ) = (0,0) and a degenerate solution 
(w s ,v s ) — (0, arbitrary) that exists for a = 0. How- 
ever, the degenerate solution does not correspond to any 
dark soliton solution, even though it is sometimes stable, 
as there is no nonlinearity in the equation for v s when 
w s = 0. 

To relate the solutions ([?]) to a specific physical situa- 
tion requires that the plane-wave variables r and 77, and 
hence f2 and A, be known, since r specifies the solutions 
and 77 is scaled into the normalized units and grating pa- 
rameters. Now f2 is determined by the angle of the plane 
wave to the Z-axis and is zero at normal incidence. In 
contrast, A is not directly accessible in experiments; it 
depends on the intensities of the fundamental and SH 
and must be determined through a nonlinear dispersion 
relation. In the limit 7 < 1, the cubic terms can be 



ignored, which gives r = (fi' ± \J p' 2 + 12/|p| 2 )/6, where 
-T=|u)| 2 + |w| 2 is the averaged intensity. In unsealed units 

r\r = (/3± a//? 2 + 12/|pun| 2 )/6, where pun = VP has now 
no dependence on 77. Which root of the radical is taken 
depends on the relative size of the two field intensities. 
When cubic terms need to be included, a simple analytic 
expression for r or rr\ is not possible and the calculation 
must be done numerically. As an example, Fig. ^ illus- 
trates the dependence of rr\ on / and (3 for a symmetric 
nonlinear grating. The two branches of solution are dis- 
tinguished by different ratios of the fundamental to SH. 
Where a choice of 77 is required in Sees. IV -VI, we sim- 
ply take 77 = 1, which can be satisfied by an appropriate 
choice of intensity / and input width Xq. 



Following standard MI analysis |T^,^,Q, we assume 
that the plane-wave solutions w s and v s of Eq. (^|) are 
perturbed by small modulations: 

w(z) — w s + 61 (z) exp (-ira) + 82 (z) exp (ivx) 

v(z) = v s + (^3(2;) exp {—ivx) + 5l(z) exp (ivx), (9) 

where * denotes the complex conjugate. The linear 
growth is governed by d z A(z) = MA(z), where 
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a bed 
—b —a — d — c 
2c 2d g 
-2d -2c -g 



(10) 



with a = -v 2 /2 — r + 7(|w s | 2 — 2w 2 s ), b = pv s — ^w 2 , 
c = pw s + ~fv s w s , d — jVgWs and g = — v /A — a + 277ii 2 . 
For MI to be absent, all the eigenvalues of M must be 
imaginary. In practice, we calculate the gain at a finite 
number of values of v, out to a maximum of v = 50. 

For r = 1 and r — 0, stability analysis shows that all 
nondegenerate and nontrivial solutions are unstable. For 
the r = — 1 solutions, however, the (LI) branch of so- 
lutions is stable. Figure |] shows the region of stability, 
which differs significantly from that found in Ref. (!§] 
because of the asymmetric nature of the induced cubic 
nonlinearity in Eq. (^|) . The continuous line in Fig. ^ is 
the boundary of (LI). Now on this branch, 7 and a can- 
not both be small, suggesting that there are two possible 
ways to achieve stability. One is to have a relatively large 
cubic nonlinearity 7, and the other a large mismatch a. 
Increasing the size of 7 basically entails decreasing p, 
the size of the effective quadratic nonlinearity, if the per- 
turbation theory is to remain valid. However, in most 
applications a large effective is desirable. The other 
possible route to stability is to increase a, i.e. to increase 
— /3' — —j3/r). However, this large residual mismatch may 
also invalidate the perturbation theory. (77 cannot be re- 
duced, because this would decrease 7, thereby moving the 
solution further from the stable region.) Thus whether 
this region of stability truly exists must be checked by 
more exact methods. 
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FIG. 3. Parameter 7-77 versus averaged intensity I for 
P = -10 (solid), (3 = (dashed), /3 = 10 (dotted), with 
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FIG. 4. (a) Region of stability for the r = 



~ 1 solutions. 



FIG. 6. Gain profiles of the r = 



— 1 solutions for 7 = 0. 



Figure plots the maximum gain of the unstable 
branches (G,U1,U2). It shows several regions where the 
gain is small, namely close to the existence boundaries 
of the solutions (G,U1,U2) and also for (G) when 7^1, 
a <C —1. Most of these 'almost stable' regions corre- 
spond to either w s or v s becoming small. For example, 
near a = 0, the solutions (G) and (U2) approach the de- 
generate solution (w s = 0) for 7 < 1/4; and when 7> 1, 
a <C — 1, the SH in (G) approaches zero. In these cases, 
the averaged equations effectively become either linear 
or a NLS equation. The exception is the branch (Ul): 
its gain becomes small in the region near where it meets 
the stable (LI) branch, as expected. Now some of these 
solutions have such a low gain that for experimental pur- 
poses they will appear to be stable. Figure || tells us that 
such experimentally stable solutions lie in diverse regions 
in parameter space and thus may be more accessible than 
the mathematically stable solutions. 
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(b) 
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FIG. 5. Maximum gain of the unstable branches (a) G and 
(b) Ul and U2 of the r = — 1 solutions (Black is largest gain). 
The hatching indicates the absence of a solution. 

We show in Fig. ^ the gain profiles for 7 = of the 
(G) solution, which corresponds to when the induced cu- 
bic terms are ignored, i.e. when the first-order terms 
are not taken into account in the perturbation theory. 
The same gain curves appear when there is no grating, 
just with the effective mismatch (3 replaced by the (very 
much larger) inherent mismatch [3 in the expression for 
a, which gives a larger \a\. Thus to zeroth order, the 
effect of the periodic grating is to increase and broaden 
the gain, because a much smaller \a\ is achieved than in 
the non-phase-matched case without a grating. However, 
the more accurate first-order theory involves a further ad- 
justment to the spectrum, with the inclusion of the cubic 
terms (so that 7 ^ 0). The cubic terms may decrease 
the maximum gain, perhaps even to zero for particular 
gratings. 



C. Cubic solutions: p — 

For some combinations of linear and nonlinear grat- 
ings, the cubic nonlinearity may dominate the quadratic. 
In the extreme case, the effective quadratic nonlinearity is 
zero (p = 0), caused by the linear grating (multiplied by 
the strength of the DC nonlinearity) cancelling the con- 
tribution of the nonlinear grating, which can occur for 
realistic physical parameters. The averaged-field equa- 
tions reduce to two NLS equations, which support the 
stable dark solitons (for r = —1): 

w = tanli(x)/V7, v = (7 > 0) (11) 
w = tanh(x)/<s/\4r/\, v = ± V5/|4 7 |tanh(x) (7 < 0), (12) 

where a = —1/2 in the second solution (|l2]). In the ab- 
sence of the parametric term, the two averaged fields are 
coupled only by phase-independent XPM terms, which 
is why these solutions do not possess the 'natural' \^ 
phase relation: in the first solution (|ll|), a is arbitrary, 
and the second solution ( |l2| ) only exists for a particular, 
nonzero value of the residual mismatch /?' = —3/2. 

The plane-wave backgrounds of the solitons (|ll| ) and 
( |l2| ) form part of a larger set of plane- wave solutions that 
exist in the absence of the quadratic terms: (w s ,v s ) = 
(V-r/7,0) and {w s ,v s ) = (^727^^ + ^/27). 
The v s = solution is stable for 7 > (r = —1) and 
unstable for 7 < (r = 1), whereas the v s ^ solution 
is always unstable. 

Figure [7(a) shows a simulation of the original field 
equations (2) with initial condition corresponding to the 
stable NLS soliton (|TT|) . No visible shedding of radia- 
tion occurs, and the fundamental intensity contains only 
rapid, small oscillations superimposed on the large av- 
erage beam, which is typical for stable, bright QPM 
solitons p3| , p0| ]. The soliton propagates stably beyond 
z = 150, which corresponds to 150 diffraction lengths or 
over 2000 coherence lengths. Thus the soliton is stable 
according to any reasonable experimental definition. In 
contrast, simulations of the two- wave dark soliton Jl2] ) 
over short distances (Fig. 0(b)) confirm that it is indeed 
modulationally unstable. We note that a small modula- 
tion (of spatial frequency v ~ 10) finally develops in Fig. 
|7 (a), which indicates a weak instability not picked up by 
the averaged-field equations. 
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(b) 



FIG. 7. The evolution of the fundamental intensity w, cor- 
responding to the solutions (a) Eq. ( |ll| ) with 7 = 0.08 and 
(b) Eq. (113) with 7 = —0.08. The gratings parameters are 
d'o/d? = 2.8 and a'/n = 0.36. r = -1. 



IV. EXACT PLANE- WAVE SOLUTIONS 



The stability analysis with the averaged-field equations 
found regions of stability. But these regions possibly 
lie at or beyond the limits of validity of the first-order 
perturbation theory. Additionally, propagation simula- 
tions showed weak instabilities not predicted by the av- 
eraged theory. Clearly, a thorough investigation requires 
a more exact technique that does not rely on the pertur- 
bation theory. Therefore, we numerically find the exact 
plane- wave solutions and apply an exact Floquet analysis 

f, fl5| to determine their stability. We begin with Eqs. 
, substituting in the Fourier expansions (||) and (jj). 
This gives a set of coupled equations for the plane-wave 
coefficients w n and v n : 

- (uk + r)w n + ^ D n+l-v-l w l v P = °> 
l,p 

-(raw + a)v n + a'y]g n ^ivi + 2J D n -i- p+ iwiw p = 0, (13) 
l i, P 
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FIG. 8. (a) Ratio R = \v \ 2 /\ wo\ versus mismatch pa- 
rameter a calculated exactly (asterisks) and by the aver- 
aged-field theory (solid), (b) Fourier components of solution 
for k — 700, a = 1 calculated exactly (solid) and by the aver- 
aged-field theory (asterisks for w n and circles for v n ). r = — 1. 

Figure ||(b) gives the Fourier coefficients |to n | and |u n | 
for a representative solution. Even though the ampli- 
tudes do not decrease exponentially with n, the decrease 
is sufficiently rapid such that the tenth component is al- 
most three orders of magnitude less than the DC value. 
Thus in our calculations, we typically include up to the 
components n — ±16 components in the Fourier expan- 
sions. 

That the approximate theory accurately predicts the 
averaged properties of the system is well known. How- 
ever Fig. ||(b) shows that this approach does not properly 
account for the higher-order Fourier components (partic- 
ularly the odd components of w and even components of 
v), which are important for a correct MI analysis 



where D n =d' g n +d' S nj o. We may rewrite these equations 
more succinctly as F„ = 0. 

We find the solution to Eq. (|1|) numerically, using a re- 
laxation technique based on Newton's method. An initial 
guess of the components w n v n (which we represent by 
the vector u^), calculated from the approximate first- 
order solutions, is iteratively corrected by the increment 
Au( fc+1 ) = - J _1 (u( fe ))F(u( fe )), where the Jacobean ma- 
trix is 



{Jh, q = 



dFi OF, OF, dFi 

8w q dv q dw* dv* 

OF* OF* dF* dF t ' 

8w„ 8v„ 8A* dB* 



(14) 



Figure g shows some of the exact solutions for r = — 1 
when there is both a linear and nonlinear grating present, 
with parameters: a'/n = 0.158, d' /d' — 5/3 and 77 = 1. 
These gratings correspond to the GaAs/GaAlAs struc- 
ture reported in Ref. p9j . Figure ||(a) shows the ratio 
^ = l u o| 2 /|wo| 2 versus the mismatch parameter for 
grating wave numbers k — 700 and K — 100. These val- 
ues of k give 7 < 0.1, which means that these solutions 
lie outside the stability region predicted by the average 
theory (Fig. ||) . Not shown are the solutions at a = for 
which w is zero or negligibly small. As the figures show, 
in terms of the DC intensity ratio, the exact solutions 
consistently agree well with the approximate solutions 
found by averaging. 



V. EXACT STABILITY ANALYSIS 

A. Floquet theory 

We determine the linear stability of the plane-wave so- 
lutions (w s ,v s ) by using exact Floquet theory |lf|. For 
perturbations of transverse frequency offset v: 

w(x, z) — w s (z) + Si(z) exp (— ivx) + S^iz) exp (ivx), 
v(x, z) = v s {z) + #3(2) exp {—ivx) + 5^(z) exp {ivx), (15) 

the linear evolution is governed by the linearized equation 
d z P = M{z)P, where 









*- 




, M = i | 









a b c* 

-b* -a -c 

2c d 

-2c* -d 



(16) 



with a = —v 2 /2 — r, b — v s {z)x{z) exp {inz), c = 
w s {z)x{z) exp (— ikz), and d = —v 2 /4 — a + a! {z). Be- 
cause of the periodicity in the coefficients a' and \' an d 
in the solutions w s and v s , the matrix M governing the 
growth of linear perturbations is also periodic, with pe- 
riod z p = 2it/\k\. To determine the stability, we need 
only find the eigenvalues of the mapping ST of the 
solution over one period P(z + z p ) = STP{z). 

Now for {P(z)}i = Si,k for some k S {1, 2, 3, 4}, the so- 
lution one period later P{z + z p ) will be the &th column 
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of ST. Thus we may construct the constant matrix ST by 
integrating the linearized equation <9 Z P = M(z)P out to 



for the four initial conditions {P(zo)}i 



corresponding to each k S {1,2,3,4}. This is done nu- 
merically with a fourth-order Rungc-Kutta method. If 
all the Xi lie on the unit circle for every v, then the solu- 
tions w s and v s are stable; otherwise there is a net gain 
whose profile is well-estimated by the maximum eigen- 
value growth rate g — max{5R[ln(Ai)/z p ]}. 



B. Applications I: Efficient phase matching 

To help in the analysis of the Floquet results, we divide 
the situations into two regimes: efficient phase matching, 
where the residual effective mismatch is much smaller 
than the material mismatch (/?' <C /?')? an d poor phase 
matching. In the former case, the averaged-field theory 
is expected to hold, whereas in the latter the large resid- 
ual mismatch invalidates the assumptions made in the 
perturbation theory (i.e. j3' <C k). 

We consider first the case which corresponds to having 
no effective quadratic nonlinearity in the averaged-field 
analysis (i.e. p = 0), to see whether the Floquet anal- 
ysis confirms the strong prediction of the averaged-field 
theory that the solutions and their MI properties are gov- 
erned by cubic nonlincarities. 
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FIG. 9. Gain curves for the dual plane-wave solution when 
(a) p = and (b) \p\ = 2/(15tt), for 7 < and a = -0.5 
(f3' — —3/2). The dashed curve gives the averaged-field re- 
sults, the solid curves give the exact Floquet results. 

To make the comparison with the exact Floquet anal- 
ysis, consider a physical system with the parameters: 
a'/n = 5/14 and d' /d' = 14/5, for which p — 0. Figure 
^(a) shows the MI gain profiles (gain versus the spatial 
frequency v of the perturbation) for the v s 7^ solution, 
corresponding to the unstable background in Fig. 0(b), 
calculated by both the exact Floquet theory (solid line) 
and the averaged-field theory (dashed) . The Floquet the- 
ory predicts that each of the two solutions has its own 
distinct gain curve, whereas in the approximate calcu- 
lation, the two solutions (±w) give the same degenerate 
gain curve. The approximate and exact results agree 
very well away from v — 0, but noticeably differ near 
v = 0. For comparison, Figure ^|(b) plots the gain curve 
for \p\ = 2/(157r), using a! / k = 1/3. Here the Floquet 
curves for the two solutions differ more markedly, but 
each agrees well with its approximate counterpart, which 



now correspond to points on the (U2) and (G) branches 
in Fig. ||. The discrepancy in Fig. g(a)in is due to higher- 
order terms neglected in the averaged theory which break 
the degeneracy of the NLS solution (|l2|). These terms, of 
order e 2 and higher, would need to be included for con- 
sistent theory that accurately predicts the gain when v is 
smaller than order e , and especially when the quadratic 
(order e ) nonlinearity is absent. 

The v s = solution is predicted by the averaged-field 
theory to be stable. However the exact Floquet theory 
shows that there is some nonzero gain at high spatial 
frequencies. However, the gain curves are extremely nar- 
row and have low peak values: the largest region of gain 
at v = 9.9 is only Av = 0.05 wide with a maximum of 
g = 0.21, which corresponds with the weak instability of 
the simulation in Fig. 7](a). 

A Floquet analysis of the solutions shown in Fig. || does 
not find any points that are stable. This is not surpris- 
ing, because these points lie outside the stability region 
predicted by the first-order analysis. But exact solutions 
corresponding to the stable lowest branch of averaged so- 
lutions can be found for lower values of the grating 
wave number. A plot of the maximum MI gain ([iT)|(a)), 
however, shows that these solutions exhibit a weak, but 
nonzero MI gain. This residual gain (~ 0.2 — > 1) in- 
creases for decreasing n and is caused by narrow gain 
bands that are once again present at high v. Other sim- 
ulations also confirm that such high-^ peaks, which are 
not predicted by the averaged theory, are a general fea- 
ture of the gain spectrum for efficient QPM gratings. 




FIG. 10. Maximum gain for solutions with (a) 5 < re < 30 
(b) 1 < re < 10. The region of zero gain predicted from the av- 
eraged-field theory is shaded. d! /d! — 5/3 and a'/re = 0.316. 

To analyze the MI spectrum in detail and illustrate 
the physical origin of these high-u peaks, we consider the 
spectra of two concrete examples that are typical for effi- 
ciently phase-matched gratings (|^| <C |/3|). Figure |ll](a) 
is the spectrum for a conventional domain-reversal grat- 
ing in LiNb0 3 (a' = d' = 0) for r=-l, n= 100 and 
exact phase matching (a— —2), and Fig. |ll|(b) is that 
of a GaAs/AlAs superstructure with a nonlinear grating 
etched though quantum- well disordering |37j (d' G /d'=A.6) 
and a linear grating chosen to be o//k=10/46. The for- 
mer is unstable according to the averaged-field analysis 
and possesses a low-^ gain that is accurately predicted 
by the average theory (Fig. |ll](c)); the latter lies in the 
predicted stable region and thus has no gain at low v. 
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FIG. 11. Floquet gain spectra for (a) the LiNbOs and 
(b) the GaAs/AlAs gratings, for a— — 2. The diamonds 
give the averaged-field results in (c) and the equivalent 
non-phase-matched homogeneous results in (d). 

The high-z/ gain bands are related to the inherent MI in 
the non-phase-matched, homogeneous \^ material (i.e. 
with no grating), as demonstrated by the close match be- 
tween the homogeneous gain band (Fig. fy](d)) and one 
of the peaks in Fig. |ll](b). In general, each gain peak 
that appears in the homogeneous spectrum usually also 
appears in the Floquet spectrum, typically with the same 
height and spectral location. Often, as in Fig. |ll|(b), the 
'homogeneous' peak is split into several closely spaced 
secondary peaks by the grating. Sometimes, as for a sym- 
metric grating (d' =0) such as in LiNb0 3 (Fig. 0(a)), the 
original homogeneous peak is suppressed, leaving only 
the secondary peaks in the high-z/ part of the spectrum. 



a nonzero gain which is large and broadband. The ex- 
act gain curve in Figs. [l^(b) (re = 6, a = 1) consists 
of a number of narrow peaks distributed over a range 
of frequencies, whereas the first-order analysis predicts a 
larger, broadband curve that lies at low frequencies. 
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FIG. 12. Gain curve for (a) re = 2 and a — 4 and (b) 
re = 6 and a = 1, calculated from Floquet theory (solid) and 
simulations (dotted). The averaged theory predicts a stable 
solution in (a) and a gain marked by the dashed line in (b). 
d' /d' = 5/3 and a'/re = 0.316. 



In Sec. VB 



where the grating efficiently phase- 
matches the fields, we saw that the MI gain spectrum falls 
into two, well-separated parts: the low-frequency bands 
induced by the periodicity and accurately predicted by 
the averaged equations, and the high-frequency bands re- 
lated to inherent x'' 2 '' gain. When the residual mismatch 
\0\ is larger (and not so different from |/3|), we find that 
the two regions no longer are well separated and that 
the averaged equations no longer accurately predict the 
low-^ bands. For very poor phase matching as in this 
section, the two features mix, producing a complicated 
gain spectrum (e.g. Fig. |l2](b)) that cannot be simply 
related to the distinct underlying physical mechanisms. 



D. Experimental Stability 



C. Applications II: Poor phase matching 



As we saw in Sec. Ill E , attempting to reach the stable 
branch of the approximate solutions while maintaining 
a large quadratic nonlinearity (p) usually involves the 
poor-phase-matching regime, where the averaged theory 
no longer holds. To illustrate what happens, we show in 
Fig. ^(b) the maximum growth rate of instabilities for 
the grating d' /d' — 5/3 and a'/n = 0.316. Here — $' > 3 
is rather large and re < 10 is fairly small, so the first- 
order perturbation results are not guaranteed to be valid 
in this region. That the Floquet theory found solutions 
with large gain in the predicted stable region is therefore 
not surprising. 

The differences between the averaged and exact treat- 
ments for these parameters are further illustrated by 
comparison of the gain profiles, which we do in Fig. [l^ 
for two representative points. The averaged solution cor- 
responding to Fig. |l2|(a) (re = 2, a = 4) lies in the 
predicted stable region, but the Floquet analysis gives 



In all our results, we found no solutions that, mathe- 
matically speaking, were modulationally stable. Even in 
the region where the averaged-field analysis predicts sta- 
bility, there is a residual gain corresponding to the high-i' 
bands. However this residual gain is not large (<? < 0.3 or 
smaller for material mismatch \(3'\ ~ 10 or larger), and 
may be neglected under a reasonable definition of exper- 
imental stability. If perturbations do not increase by any 
orders of magnitude within a few diffractions lengths, 
which is a typical experimental length scale, then the in- 
stability will not be detectable. However if we do neglect 
these high-z/ bands, so that the averaged and Floquet re- 
sults coincide, then we must also allow some solutions 
previously classed as unstable by the averaged-field anal- 
ysis to be classed as stable, because the \ow-v bands are 
sometimes smaller than the high-i^ bands. These addi- 
tional solutions correspond to the light-grey regions in 
Fig. | for r = — 1 as well as similar regions for r = and 
r = 1 , which were previously regarded as totally unstable 
branches. 
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VI. PROPAGATIVE SIMULATIONS 

To confirm the Floquet predictions of the MI spec- 
tra, we perform numerical simulations of the (scaled) 
field equations (Q). We can determine gain curves by 
launching plane waves seeded with small fluctuations, 
which then evolve linearly. Using Gaussian white noise 
to excite all frequencies, we calculate the gain curve 
from the Fourier transformed amplitude J3(| E± as: 
g(y) = [ln\E l (Z 2 ,T)\-ln\E 1 (Z 1 ,T)\}/(Z 2 -Z 1 ), where 

T = y/fjV. 

Propagation calculations confirm, for example, the 
gain curves predicted by Floquet theory in Fig. g(a). Im- 
portantly, they agree with the Floquet curves in regions 
where the averaged-field and Floquet theories predict dif- 
ferent results, such as near v = in these figures. A more 
dramatic illustration of this is shown in Figs, [l^, where 
the averaged and exact predictions markedly disagree. 
Small, narrow peaks are difficult to distinguish from the 
noise floor in these simulations, but the peaks seen in the 
smoothed data in Fig. |l2|(b) (dotted) clearly correspond 
to the largest Floquet peaks (solid). For the same rea- 
son, the small, narrow high-^ bands predicted by Floquet 
theory are difficult to detect, but with a careful choice of 
simulation parameters these peaks can be seen, as shown 
in Fig. for the LiNb03 grating. Clearly, the simula- 
tions confirm the Floquet calculations, in both the low-v 
[Fig. |l^(b)] and high-^ [Fig. |l3|(c)] regions. 
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FIG. 13. (a) Gain profile calculated from propagative sim- 
ulations, for the same parameters as in Fig. |ll](a). Compar- 
isons with the Floquet theory (solid curve) are given in (b) 
and (c). 



VII. CONCLUSIONS 

In summary, we determined the exact structure and 
modulational stability of plane-wave solutions in peri- 
odic x^ 2 ) materials, using both Floquet transfer-matrix 
calculations and direct numerical simulations. In partic- 
ular, we found Fourier expansions of the plane waves in 
the propagation direction and used these in calculations 
of the gain spectrum (growth rate versus transverse fre- 
quency of the instability). In the regime of efficient phase 



matching, we found that, due to the periodicity, the spec- 
tra contain gain bands both at low spatial frequency, 
which are caused by the phase-matching gratings, and 
at high spatial frequency, which are related to the inher- 
ent gain of homogeneous materials. For poor phase- 
matching, the gain bands do not fall into these different 
regions of the spectrum and so cannot be respectively 
attributed to distinct underlying causes. 

In the regime of efficient phase matching, the exact 
calculations were supplemented by a simpler averaged- 
field approach characterized by effective quadratic and 
induced cubic nonlincaritics. For sufficiently large cu- 
bic nonlinearities and away from exact phase match- 
ing, the averaged-field theory predicts stability for some 
solutions. However, the exact techniques showed that 
whereas the averaged-field approach accurately predicts 
slowly varying properties of the system, it sometimes 
inaccurately treats the higher Fourier components, as 
these are treated to a lower order in the theory. It 
also shows inaccuracies in the MI predictions, being un- 
able to predict the high-frequency gain bands. Neverthe- 
less, the averaged-field theory accurately describes low- 
frequency gain, even in the extreme case where the ef- 
fective quadratic nonlinearity is cancelled by competi- 
tion between the linear and nonlinear gratings. Futhcr- 
more, because the high-frequency gain is often small, 
the averaged-theory zero-gain predictions can be accu- 
rate under experimentally relevant definitions of stabil- 
ity. 
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